// Interpolator.cpp - MARIN MSG Anneke Sicherer-Roetman 20081105

#include <math.h>    // cos
#include <algorithm> // std::max and std::min
#include "Interpolator.h"
#include "assertplus.h"

#define PI 3.14159265359

/// constructor
/** Creates table with Chebyshev nodes. */
Interpolator::Interpolator(const Array1D<double> &x, ///< original x coordinates
                           const Array1D<double> &y) ///< original y coordinates
    : m_x(x)
    , m_y(y)
    , m_t(x.size())
{
    ASSERT(x.size() == y.size(), "");
    int nn = x.size();
    double xmin = 1e300;
    double xmax = 0;
    for (int i = 0; i < nn; ++i)
    {
        xmin = std::min(x[i], xmin);
        xmax = std::max(x[i], xmax);
    }
    for (int i = 0; i < nn; ++i)
    {
        m_t[i] = 0.5 * (xmax + xmin) + 0.5 * (xmax - xmin) * cos (PI * (2 * i + 1) / (2 * (nn - 1) + 2));
    }
}
   

/// destructor
/** Does nothing. */
Interpolator::~Interpolator()
{
}


/// performs Chebyshev interpolation
/** The interpolation is done for the desired number of points. */
void Interpolator::interpolate(int n,               ///< number of points to interpolate
                               Array1D<double> *px, ///< interpolated x coordinates
                               Array1D<double> *py) ///< interpolated y coordinates
{
    px->resize(n);
    py->resize(n);
    int nn = m_t.size();
    double step = (m_t[nn - 1] - m_t[0]) / n;
    double point = m_t[0];
    for (int tp = 0; tp < n; ++tp)
    {
        double xx = 0;
        double yy = 0;
        for (int k = 0; k < nn; ++k)
        {
            double weight = 1;
            for (int j = 0; j < nn; ++j)
            {
                if (k != j)
                {
                    weight *= (point - m_t[j]) / (m_t[k] - m_t[j]);
                }
            }
            xx += m_x[k] * weight;
            yy += m_y[k] * weight;
        }
        (*px)[tp] = xx;
        (*py)[tp] = yy;
        point += step;
    }
}

